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We numerically simulate mechanically stable packings of soft-core, frictionless, bidisperse disks in 
two dimensions, above the jamming packing fraction (fj. For configurations with a fixed isotropic 
global stress tensor, we investigate the fluctuations of the local packing fraction <))(r) to test whether 
such configurations display the hyperuniformity that has been claimed to exist exactly at <j)j. For 
our configurations, generated by a rapid quench protocol, we find that hyperuniformity persists only 
out to a finite length scale, and that this length scale appears to remain finite as the system stress 
decreases towards zero, i.e. towards the jamming transition. Our result suggests that the presence of 
hyperuniformity at jamming may be sensitive to the specific protocol used to construct the jammed 
configurations. 
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I. INTRODUCTION 

When a system of athermal (T = 0) particles with 
only contact interactions is compressed, it seizes up into 
a rigid disordered solid at a critical value of the packing 
fraction (pj known as the jamming transition [T]-^. For 
a system of monodisperse frictionless spheres at pj, it 
was observed numerically HE] that density fluctuations 
appear to be suppressed on long length scales, with a 
structure function 5'(q) (density-density correlation) that 
vanishes as S'(q) ^ |q| when the wavevector |q| —^ 0. 
This is in contrast to behavior in a normal liquid where 
S'(q —>■ 0) —?► constant. Such a system with suppressed 
density fluctuations has been denoted as “hyperuniform” 
0 - 

However when spheres that were bidisperse or polydis- 
perse in size where studied, this characteristic feature of 
S'(q) was no longer observed, and ^(q —>■ 0) was found 
to be finite H H. It was then argued by Berthier et 
al. and by Zachary et al. uniEi, that in such size- 
disperse systems it is the fluctuations of the packing frac¬ 
tion (p, rather than fluctuations of particle density, that 
are suppressed at pj (for monodisperse systems, packing 
fraction fluctuations and density fluctuations become the 
same at long wavelengths). The presence of such hyper¬ 
uniformity of the packing fraction at jamming would be 
important, as it would provide a purely structural means 
for distinguishing particles in a disordered jammed con¬ 
figuration from those in a liquid, and perhaps provide a 
way to determine a diverging length scale as the jamming 
transition is approached [T^ . 

In this work we consider mechanically stable packings 
of bidisperse, soft-core, frictionless, disks in two dimen¬ 
sions at finite isotropic global stress above the jamming 
transition pj. Our configurations are generated by a 
rapid quench protocol. We test these configurations for 
hyperuniformity using both real-space and wavevector- 
space methods. We find that hyperuniformity persists 


only out to a finite length scale, and that this length 
scale appears to remain finite as the system stress de¬ 
creases towards zero, i.e. as one approaches the jamming 
transition. Moreover, we argue that measuring fluctua¬ 
tions at a given wavevector q gives a better test of hyper¬ 
uniformity than measuring fluctuations over a real-space 
window of length R, as the latter can be strongly effected 
by the fluctuations on all length scales smaller than R, 
whereas the former measures fluctuations specifically on 
the length scale 2TT/q. 

The remainder of this paper is organized as follows. 
In Sec. we define what we mean by the local pack¬ 
ing fraction ^(r) and discuss the wavevector-dependent 
and real-space measures we will use to test for hyper¬ 
uniformity. In Sec. m we describe the details of our 
numerical model and the minimization method we use 
to construct mechanically stable configurations at fixed 
isotropic global stress. In Sec. |TV] we present our nu¬ 
merical results. In Sec. El we discuss our results and 
make comparisons with recent works on this topic. The 
Appendix provides further details about the accuracy of 
our numerical minimization method for constructing our 
configurations. 


II. LOCAL PACKING FRACTION 

In this section we define the quantities we will com¬ 
pute in order to test for hyper uniformity. Here we define 
quantities as appropriate to a system of two dimensional 
circular disks, so as to match our numerical simulations, 
however the generalization to higher dimension or other 
shaped particles is straightforward. 

Consider a polydisperse collection of N disks in a sys¬ 
tem of total volume V, satisfying Lees-Edwards bound¬ 
ary conditions m- Disk i has its center located at posi¬ 
tion Vi and has volume Vi (in our two dimensional system 
we will use “volume” to mean area). The local particle 
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density can then be written as. 


fluctuations in the packing fraction at wavevector q are 

n(r) = E^(’' 


given by. 


(1) 

X(q) = ^(l^'q^-q)- 

(10) 

i 

Defining the Eourier transform. 


The signature of hyperuniformity is then 


Tiq = / dPr e^^'^n{v), 

Jv 

(2) 

x(q) ~ |q| as |q| 0, 

(11) 


the structure function (density-density correlation) is, 

5'(q) = ^(nqn_q), (3) 

where here, and henceforth, (...) denotes an average over 
independently quenched configurations. For the bidis- 
perse systems we study here, we expect S (q) to approach 
a constant as |q| —)■ 0 [7HII]. 

The global packing fraction of the system is defined as, 

(4) 

i 

For the local packing fraction (j){r), two slightly different 
definitions have been proposed in the literature. Zachary 
et al. nnun] use a definition that is equivalent to, 

definition! : = ^A,(r-r,), (5) 


where the indicator function (r) is such that for a par¬ 
ticle centered at the origin, 


Ai(r) 


1, if r lies within the area of the particle 
0, otherwise 


( 6 ) 


so that fy cPr Ai(ri) = vt. 

Berthier et al. [5] use a definition m that is equivalent 
to. 


definition!! : ^ Vi6{r — r^). (7) 

i 

Both definitions give correctly the global packing fraction 
ofEq. 0, 

y = ( 8 ) 

Definition ! spreads the weight of each particle uniformly 
over its area, while definition !! treats each particle as a 
point object with weight equal to its area. Definition !! 
views the particle positions as a point process, while def¬ 
inition ! views the particles as defining a heterogeneous 

medium Ha- 

Defining the Fourier transform, 

(l)q = (9) 


whereas x(q 0) —>■ constant if the system is not hype¬ 
runiform. 

Note, using definition !! of Eq. 0 we have, 

= ( 12 ) 

i 

whereas using definition ! of Eq. 0 we have, 

<(.q = ^A,qe*‘i'-' (!3) 


where A^q is the Eourier transform of Ai(r). Since A^q —>■ 
Vi as |q| —>■ 0, the two definitions of Eq. 0 and 0 must 
give the same x(q) in the limit |q| —>■ 0, hence both are 
in principle good measures for hyperuniformity. 

Note, for circular disks in two dimensions, A^q depends 
only on the magnitude |q| and is given by, 

\q = vJ{\q\di/2), f{y) = ^ f dxxJoix). (!4) 

y Jo 

Here di is the diameter of the particle i and Jo)^;) is the 
Bessel function of the first kind. 

We will also consider another wavevector dependent 
measure of hyperuniformity, as introduced by Berthier 
et al. [5], the thermal compressibility XT(q) defined by. 


definition!!!: [nTy^lq)] ^ = y^^XsS^J, {q)Xs'■ 

s,s' 

(15) 

Here n = N/V is the particle density, s and s' label 
distinct species of particles of given diameter ds, Xg = 
Ng/N is the global concentration of species s, and S~^ (q) 
is the inverse of the matrix, 

'S'ss'(q) = 'J^i'ngqngi — q) , (16) 


where rtgq is the Eourier transform of the particle den¬ 
sity of species s alone. The quantity XT(q) in Eq- (15) 
is derived as the compressibility of a polydisperse liq¬ 
uid of particles in thermal equilibrium at temperature T. 
Eor our nonequilibrium athermal system, in which fluc¬ 
tuations from configuration to configuration are induced 
by our rapid quench protocol rather than a finite tem¬ 
perature, the physical interpretation of XT(q) as a com¬ 
pressibility is unclear; nevertheless the right hand side of 
Eq. (151 is an interesting measure of density fluctuations, 
and so we will compute it for the sake of comparison. 
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We will also consider hyperuniformity as measured in 
real-space by computing the fluctuations of (/)(r) over a 
circular window of radius R. Place a circle of radius R 
at a random position within the system and denote this 
region as the volume Vr- We can then define the average 
packing fraction on this region as, 

<i^R= (17) 

The difference in (j)R between using definition I of Eq. (§ 
and definition II of Eq. 0 for (^(r) is then as illustrated 
by the sketch in Fig. In definition I we count all over¬ 
lapping volume between particles and the volume Vr] 
particles which are not entirely contained within Vr con¬ 
tribute only the overlapping fraction of their volume, as 
illustrated. In definition II we count the entire volume of 
particles whose centers lie within the volume Vr] parti¬ 
cles whose centers lie outside Vr contribute nothing, even 
if they overlap Vr. 




FIG. 1. Shaded volume represents the volume that con¬ 
tributes to (fiR according to (a) the definition I of ^(r) in 
Eq. (§, and (b) the definition II of (/)(r) in Eq. 0 . 


We then compute the variance, 

var((/)i^) = - {cj)R)^. (18) 


One can show that va,T{(j)R) is related to x(q) by 

E x(q)A"flq, (19) 


where A^q is the Fourier transform of the indicator func¬ 
tion Aij(r) for a circular volume of radius R, and the sum 
is over all q consistent with Lees-Edwards boundary con¬ 
ditions excluding q = 0 m- 

When x(q) —>• constant as |q| —>■ 0, as in a liquid, the 
above gives [101 E] for tbe limiting large R behavior in 
two dimensions. 


foiTiquid : var(0i^) - —. (20) 

For a hyperuniform system, with x(q) |q| as |q| —>■ 0, 
the limiting large R behavior in two dimensions is [lOllllj . 


a + bln R 
^3 


( 21 ) 


Since the |q| —>■ 0 limiting behavior of x(q) must be 
the same for definitions I and II, we expect that the large 
R limiting behavior of var(0^) must in principle also be 
the same. However, unlike what we will find for x(q): 
will find that for the system sizes and length scales we 
can simulate, var((/)ij) vs R behaves very differently for 
the two definitions of ((>(r). 

The relative merits of the wavevector-dependent 
method x(q) compared to the real-space method 
var((j)R), for detecting hyperuniformity as applied to par¬ 
ticle images from physical experiments, has recently been 
discussed in Ref. m 

III. MODEL 

Our two dimensional system of N particles is a bidis- 
perse mixture of equal numbers of big and small circular, 
frictionless, disks with diameters db and dg in the ratio 
db/dg = 1.4 |2]. Disks i and j interact only when they 
overlap, in which case they repel with a soft-core inter¬ 
action potential. 




^/Cg(l ‘^ij / dij^ , Tij < dij 
0 , ’^ij dij . 


( 22 ) 


Here Vij is the center-to-center distance between the par¬ 
ticles, and dij = {di + dj)/2 is the sum of their radii. 
We will measure energy in units such that fcg = 1, and 
length in units so that the small disk diameter dg = 1. 
Unless otherwise stated, our results are for the harmonic 
interaction with a = 2. 

The geometry of our system box is characterized by 
three parameters, Lx,Ly,j, as illustrated in Fig. 
and Ly are the lengths of the box in the x and y di¬ 
rections, while 7 is the skew ratio of the box. We use 
Lees-Edwards boundary conditions m to periodically 
repeat this box throughout all space. 


yLy 



FIG. 2. Geometry of our system box. Lj, and Ly are the 
lengths in the x and y directions, and 7 is the skew ratio. 
Lees-Edwards boundary conditions are used. 


In this work we consider only packings with an 
isotropic total stress tensor 

Ea/3 = rjv<5a,s, where Tn=pV, 


for hyperuniform : vai^cfiR) 


(23) 
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p is the system pressure, and V = L^Ly is the total 
system volume. Here a, [3 denote the spatial coordinate 
directions x, y. 

To construct such isotropic packings, in which the 
shear stress vanishes, we use a scheme in which we vary 
the box parameters Lx,Ly and 7 as we search for me¬ 
chanically stable states [15] • We introduce m a modified 
energy function U that depends on the particle positions 
{rj, as well as Lx,Ly,j, 


U = U + TN{lnL,+lnLy), C/ = ^V,,(r„). (24) 

i<j 


Noting that the interaction energy U depends implic¬ 
itly on the box parameters Lx,Ly,"/ via the boundary 
conditions, we get the relations. 


dU 


Lrr. ^xx “f 

dU 


dU 


= -s. 


^XX : ,^xy, ^Xy, 


(25) 


— ^yy 




Starting from an initial configuration of randomly po¬ 
sitioned particles in a square box [L^ = Ly,^ = 0 ) at 
packing fraction t/iinit = 0.84, and fixing a target value 
of Tjv, we then minimize U with respect to both particle 
positions and box parameters. Our minimization can be 
considered as a rapid quench from infinite to zero tem¬ 
perature, keeping the final total system stress fixed. The 
resulting local minimum of U gives a mechanically stable 
configuration with force balance on each particle and a 
total stress tensor that satisfies 


'^XX — Syy — TjV, '^xy — 0 . (26) 

For minimization we use the Polak-Ribiere conjugate 
gradient algorithm |20| . We consider the minimiza¬ 
tion converged when we satisfy the condition (Ui — 
Ui+ 5 o)/Ui +50 < e = 10 “^°, where Ui is the value at the 
ith step of the minimization. Tests that our procedure 
gives well minimized configurations are discussed in the 
Appendix. Our results at each value of Tat are averaged 
over 1000-10000 (depending on the system size) indepen¬ 
dently generated isotropic configurations. Configurations 
are generated independently at each value of Tat- 


IV. RESULTS 

We simulate our system for a range of total system 
stresses Tat spanning just over two orders of magnitude. 
It will be convenient to parameterize our configurations 
by the intensive quantity p = Tm/N, the total stress 
per particle; p is related to the ordinary pressure p by 
p = p{V/N). We have considered four different system 
sizes, N = 8192, 16384, 32768 and 65536, each for equal 
values of p = 0.0001373 to p = 0.0183105. We use large 
systems in two dimensions so as to be able to probe small 
wavevectors q, and so test for hyperuniformity on long 
length scales. 


A. Global quantities 


Before considering the behavior of local packing frac¬ 
tion fluctuations, we first consider several global proper¬ 
ties of the system in order to establish where our sys¬ 
tems lie with respect to the jamming transition. For our 
model, a detailed finite-size-scaling analysis | 2 T| found 
that a rapid quench from random positions at fixed pack¬ 
ing fraction gave a jamming fraction of fij = 0.84177. 
However, since it is established [551 - 125] that the jam¬ 
ming fraction fij of mechanically stable configurations 
can depend on the specific protocol used to produce those 
configurations, there is no guarantee that fij for rapid 
quenching to constant stress p necessarily results in the 
same exact value of fij. 


Since our minimization procedure varies the box 
lengths Lx and Ly to achieve the desired global stress 
Fat, different configurations at a common fixed Fat may 
have slightly different volumes (see Appendix), and hence 
different global packing fractions. In Fig. [^ we plot 
the average global packing fraction {fi) vs the stress per 
particle p, for systems with N = 8192 to 65536 parti¬ 
cles. Panel (a) shows the results on a linear-linear scale, 
where it appears that finite size effects are negligible. 
It has been predicted [5] that, for our harmonic inter¬ 
action of Eq. (22 1 , pressure scales linearly with pack¬ 
ing fraction. In our data, however, we see a small but 
clear curvature at larger p. We thus fit (solid lines in 
Fig.§ our results to {(j>) = fij oip -I- a 2 P^, regard¬ 
ing the quadratic term as a correction to scaling. This fit 
gives fij = 0.84159 ±0.00002, with the error representing 
the variation in values obtained for the different system 
sizes. 

However, to examine more closely the points at the 
smallest p, in Fig. If we plot {(p) — 0.8415 vs p on a log- 
log scale ED- We now see a definite finite size effect in 
the results for the two smallest p values, and that these lie 
noticeably below the fitted quadratic curve (solid line). 
Our above estimate of pj should therefore be taken with 
some caution. A more accurate determination of cpj, as 
well as the power-law dependence betweenp and {fi—pj), 
should take into account these finite size effects. Such an 
analysis is outside the scope of the present work. 

The jamming transition of frictionless particles is well 
determined by the isostatic condition m, where the 
number of constraints on the particles exactly equals the 
number of degrees of freedom. For frictionless spherical 
particles this condition requires that the average number 
of contacts (z) for a given particle is equal to twice the di¬ 
mensionality of the system; for two dimensions, ziso = 4. 
Numerically, this condition is found to hold quite pre¬ 
cisely provided one first excludes from the system “rat¬ 
tler” particles [5]. A rattler is any particle which is not 
at a strict local energy minimum, but may move with¬ 
out cost in energy in one or more directions. To locate 
the rattlers in our two dimensional system we loop re¬ 
cursively through all our particles removing any particle 
with less than three contacts; i.e., after an initial pass in 
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FIG. 3. (color online) Average global packing fraction {(f)) 
vs stress per particle p, for systems with N — 8192 to 65536 
particles on (a) linear-linear, and (b) log-log scales. Solid lines 
are a fit to a quadratic function {(f)) = (f)j + aip + 02 P^. 


FIG. 4. (color online) Average excess contact number {Az) = 
{z) — aiao vs stress per particle p, for systems with N = 8192 
to 65536 particles on (a) linear-linear, and (b) log-log scales. 
Solid lines are a fit to {Az) = p^{ao + niP + + ci 2 P^), 

where we hnd x « 0.546. 


which such particles are removed, we loop again through 
the remaining particles and remove any that now have 
less than three contacts, repeating this procedure un¬ 
til no additional particles are removed. The total num¬ 
ber of removed particles is then the number of rattlers. 
Removing such rattlers and computing the resulting av¬ 
erage {z) of the remaining particles, in Fig. we plot 
{Az) = {z) — Ziso, vs p, for system sizes N = 8192 to 
65536. Panel (a) shows our results on a linear-linear 
scale, while panel (b) shows a log-log scale. In this case, 
as has been noted previously [55], finite size effects are 
truly negligible for the range of p and N considered here. 
For the harmonic interaction used here, theoretical ar¬ 
guments [53] predict (Az) ^ p^^^ close to the jamming 
transition. Since our values of p extend moderately above 
jamming, {{(f))ma.yL —4'j)/4>J ~ 0.05, we fit our data to the 
form {Az) = p^{aQ+aip+a 2 P^+a 2 P^), where the polyno¬ 
mial factor is an empirical form to account for corrections 
to scaling when not sufficiently close to jamming. We find 
the value x « 0.546 ± 0.001, with the error representing 
the variation in values obtained for the different system 
sizes. A more careful scaling analysis, going to lower 
stresses p closer to jamming, is desirable before conclud¬ 
ing the exponent is truly x > 1/2. However we may note 
that a recent reanalysis [30] of the data of Ref. [28] has 
similarly found values of a; > 1/2 in both two and three 
dimensions. 

As a final measure of the global properties of our sys¬ 
tems we consider the density of states D{uj) of the dy¬ 
namical matrix of our minimized configurations [5] 1291 
El]. Expanding the interaction energy U{{ri}) to sec¬ 
ond order in small particle displacements about the en¬ 
ergy minimized configuration defines the dynamical ma¬ 
trix. The eigenvalues A of that matrix, and correspond¬ 
ing eigenvectors, determine the response of the system 
to vanishingly small elastic perturbations. Following con¬ 
vention and assuming Newtonian equations of motion for 
the response to such perturbations, the eigenvalues A are 
related to the frequencies of the normal modes of vibra¬ 
tion w by A = In Fig.jsjwe plot the density of such fre¬ 
quencies D{u!) vs w on a linear-log scale, for a stress per 


particle ranging from p = 0.0001373 to 0.0183105. Be¬ 
cause of the numerical difficulty of computing the eigen¬ 
value spectrum for large matrices, we show results only 
for our smallest system size with N = 8192 particles; 
curves at each p are averaged over 3 independent con¬ 
figurations. We see clearfy the plateau at small w, often 
referred to as the “boson peak” [5], that shows the ex¬ 
cess of low frequency modes characteristic of a marginally 
stable solid. As p decreases, the low frequency edge of 
the plateau, uj*, moves steadily to lower values and pre¬ 
sumably vanishes as p ^ 0 [5]. Our range of stress p 
is thus clearly in the region where marginal stability is 
characterizing the structure of the packing out to ever 
increasing length scales as p decreases. 



FIG. 5. (color online) Density of frequencies of small elastic 
vibrations D{u)) vs ui. Results are shown for a system of A = 
8192 particles for a range of stress per particle p = 0.0001373 
to 0.0183105. Each curve is an average over 3 independent 
energy minimized configurations. 
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B. Wavevector-dependent fluctuations 

We now consider the fluctuations of the system at fi¬ 
nite wavevectors q. For the system geometry of Fig. 
the wavevectors allowed by the Lees-Edwards boundary 
conditions have the form, q = niihi + TO 2 b 2 , where mi 
and m 2 are integers, and the basis vectors are bi = 
( 27 r/La;)(x — yy) and b 2 = {2TTjLy)y. For simplicity 
we will look at wavevectors oriented in the y direction, 
i.e. q = mb 2 = qy, with q = 2'Km/Ly for integer m 
[32] , Because each different configuration may have a 
slightly different value of Ly, since Ly is a free variable 
determined by the targeted value of Fjv, we average data 
points at common values of m; however the variation in 
Ly over different conhgurations, while finite, is in prac¬ 
tice negligible for the large system sizes we consider here 
(see Appendix). 

In Fig. we plot the structure function S{qy), that 
measures fluctuations of particle density, vs q for a system 
of fV = 32768 particles for a range of stresses p. As 
expected, we see that S{qy) saturates to a finite constant 
as q decreases, for all p. 


creases, q* decreases slightly, but at sufficiently small p, 
the data at different p appear to be approaching a com¬ 
mon curve, with a common limiting value of q* « 0.15. 
We thus conclude that as p decreases, and one approaches 
the jamming transition, our configurations display hype¬ 
runiformity only out to a finite length scale « 27r/g* « 
42. This is the main result of this work. 
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FIG. 6 . (color online) Structure function S{qy) vs q, giv¬ 
ing the fluctuation in the density of particles at wavevector 
qy. Results are shown for a system with N = 32768 par¬ 
ticles, for several different values of the stress per particle, 
p = 0.0001373 to 0.0183105. S{qy) approaches a constant 
as 5 —i- 0. The error bars shown in the figure represent one 
standard deviation of estimated statistical error. 


In Fig. we show our results for the fluctuations of the 
local packing fraction, plotting x(gy) vs q, where we have 
used definition I of Eq. (§ for the local packing fraction 
4>{r). We show results for a system of iV = 32768 parti¬ 
cles for a range of different stresses p. We see that as q 
decreases, xiiy) decreases roughly linearly in q as was ob¬ 
served previously. However when q gets sufficiently small, 
X reaches a finite minimum at a g*, and then increases 
as g —>■ 0, rather than vanishing as expected for a hyper¬ 
uniform system. Note that the limiting value x(gy —t 0) 
is increasing as p increases. 

The variation of q* with p is quite small. As p de¬ 


FIG. 7. (color online) Fluctuation in packing fraction x(gy) 
vs g, using the definition I of Eq. § for the local packing frac¬ 
tion (j>{r). Results are shown for a system with N = 32768 
particles, for several different values of the stress per particle, 
p = 0.0001373 to 0.0183105. As q decreases, we find that 
x(gy) decreases roughly linearly, but then reaches a mini¬ 
mum below which it increases, thus implying the absence of 
hyperuniformity on large length scales. As p decreases, the 
curves of xiiv) approach a common limiting curve. The error 
bars shown in the figure represent one standard deviation of 
estimated statistical error. 


One may question whether our observed behavior of 
x(gy) at small g is not some artifact of our numerical 
procedure. In the Appendix we show a careful analysis 
that this behavior is not an artifact of an insufficiently 
converged minimization procedure. Another possibility 
might be that it is a finite size effect. In Fig. we 
therefore plot x(gy) vs g (using definition I for ^(r)) 
for several different system sizes, N = 8192 to 65536, 
for our smallest stress p = 0.0001373 and for our largest 
stress p — 0.0183105. Apart from the fact that in sys¬ 
tems with larger N we can measure down to smaller g 
(since g = 2'Km/Ly), the measured x{<iy) is found to be 
completely independent of the system size. 

Finally we consider our two other wavevector- 
dependent measures of hyperuniformity, the fluctuation 
x(q) using definition II of Eq. Q for the local packing 
fraction (/{y), and the thermal compressibility nTxT(q) 
of Eq. (15), used by Berthier et al. [S], which we denote 
as “definition III.” In Fig. we plot x(9y) ^s g for defi¬ 
nitions I, II, and III, for a system with N = 32768 parti¬ 
cles at our smallest and largest values of p. While these 
quantities all differ somewhat at the larger values of g, we 
see that definitions I and II become completely equal at 
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FIG. 8. (color online) Fluctuation in packing fraction xiiv) 
vs q, using the definition I of Eq. © for the local packing frac¬ 
tion ((»(r). Results are shown for several different system sizes, 
N = 8192 to 65536, for our smallest and largest values of the 
stress per particle, p — 0.0001373 and 0.0183105. The error 
bars shown in the figure represent one standard deviation of 
estimated statistical error. 


smaller q, in particular about the minimum q*, as should 
be expected from the discussion following Eq. (13). Def¬ 
inition III for uTxt is completely equal to definitions I 
and II at small q about the minimum q* at our lowest 
p = 0.0001373. For the largest p = 0.0183105 we find 
a small deviation between uTxt and x that persists to 
low q at and below q*; however the qualitative behav¬ 
ior remains the same. We thus conclude that all three 
approaches lead to the same conclusion: that hyperuni¬ 
formity extends only out to a finite length scale for our 
mechanically stable packings above the jamming transi¬ 
tion, and that this length remains finite as the jamming 
transition is approached. 


C. Real space fluctuations 


In this section we consider the real space fluctuations of 
the local packing fraction, defined over a circular window 
of radius R, by computing the variance of (pR as dehned 
in Eq. For each conhguration we use several dif¬ 

ferent, non-overlapping, circular windows at each given 
R. When the diameter 2R is roughly equal to half the 
length of the system E/2, we take only a single window 
per configuration. 


In Fig. 10 we plot var((/fl) vs i?, comparing results 
from using definition I of Eq. ([^ for the local packing 
fraction (j){r) with that of definition II of Eq. ([^. We 
show results for our smallest stress p = 0.00011373 and 
our largest p — 0.0183105, for a system with N = 32768 
particles. Although the corresponding x(gy) for these 
two definitions were shown in Fig. to be essentially 
identical at small q, we see a rather dramatic difference 
in the behaviors of the corresponding var((/jj) for the en¬ 


FIG. 9. (color online) Comparison of the different wavevector- 
dependent measures of hyperuniformity: “definition I” and 
“definition H” refer to the fluctuation of local packing fraction 
x(q) with b(r) determined from the definitions of Eq. (§ 
and Eq. 0 respectively, while “definition III” refers to the 
quantity nTxT(q) defined in Eq. (15). Results are shown 
for a system with N = 32768 particles, for our smallest and 
largest values of the stress per particle, p = 0.0001373 and 
0.0183105. The error bars shown in the fignre represent one 
standard deviation of estimated statistical error. 


tire range of R we study. As expected from Fig. the 
fluctuations for definition I are smaller than for definition 
II. However the two definitions also appear to give dif¬ 
ferent power-law dependencies for the decay of var((/fl) 
with R. Definition II gives roughly a 1/R^ decay, while 
definition I seems to be closer to a 1/i?^ decay at large 
length scales. We will see below that the big difference 
in magnitude of var((/)fl) comparing definition I and defi¬ 
nition II, as well as the apparent difference in power-law 
decay, can be attributed to the contributions to var((/fl’) 
from moderate to large |q| fluctuations (i.e. small length 
scale fluctuations), and that these higher |q| fluctuations 
are much larger for definition II. 


To examine this decay more c lose ly, w e c onsider 
i?^var((/i^), which according to Eqs. (20) and (21) should 
approach a constant for a liquid-like system, and (a -I- 
bhiR)/R for a hyperuniform system. In Fig. 11 we show 
R^va.i{(j>R) vs R using definition I, for our smallest and 
largest stresses, p = 0.0001373 and p = 0.0183105, for 
several different system sizes from N = 8192 to 65536. 
At small R we see that i?^var((/)fl) decays as R increases. 
A power-law ht to the small R data in panel (a) gives 
a decay ~ while in panel (b) we find ~ it 

is not clear that these exponent values have any funda¬ 
mental significance. However as R increases, this decay 
is cutoff at a length R* where i?^var(^/j) reaches a min¬ 
imum. Comparing panels (a) and (b) we see that R* 
decreases only slightly as p increases over the two orders 
of magnitude. At the lowest stress, R* « 18, correspond¬ 
ing to a window of diameter 2R* = 36; this is roughly 
consistent with the value of £* = 271 jq* « 42 obtained 











FIG. 10. (color online) Comparison of the fluctuation in local 
packing fraction, Yax{<pii) vs R, using definition I of Eq. B 
and definition II of Eq. 0 for the local packing fraction (?i>(r). 
Results are shown for a system with N = 32768 particles, 
for our smallest and largest values of the stress per particle, 
p = 0.0001373 and 0.0183105. Thick solid lines represent the 
dependancies 1/R^ and 1/R^, as indicated. 


from the minimum of x{.QY) iii Fig- 0 

For R > R* we see that i?^var(^fl) increases, rather 
than saturating to a constant as might be expected. This 
is the real space manifestation of the increase in xiQV) 
q decreases below q*. Whether R^vai{(j)fi) will continue 
to increase, or saturate to a constant, as R increases fur¬ 
ther (i.e. whether xi^v) continues to increase or satu¬ 
rates to a constant as g —> 0) remains unclear. The fi¬ 
nite size dependence seen at large R is another reflection 
of the increase in xi^f) Q decreases below q*. From 
Eq. (191 we have that i?^var(0fl) is related to the sum of 
x(q) over all allowed wavevectors. As N increases, the 
smallest allowed q decreases (gmin ~ 1/F), and we get ad¬ 
ditional contributions to this sum, resulting in the finite 
size effect at large R. That this effect is more noticeable 
at the higher stress p (compare Fig. with [TT^) is a 
consequence of the fact that the increase in x(q)^ small 
|q| becomes steeper at larger p (see Fig.[^. If x(q) even¬ 
tually saturates to a constant as |q| —>■ 0, these additional 
contributions as N increases will become a negligible part 
of the sum, and the finite-size-effect will similarly become 
negligible. 

we similarly plot i?^var((/)jj) vs R, but now 
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In Fig. 

using definition II of Eq. Q for the local packing fraction. 
Again we show results for several different system sizes, 
N = 8192 to 65536 for our smallest stress p = 0.0001373, 
and largest stress p = 0.0183105. The results are dra¬ 
matically different from what is seen in Fig. El Here we 
see a much weaker dependence on the stress p, only a 
small finite size effect at the largest R, and a clear R~^ 
decay over much of the range of data (a power-law fit 
to the data gives more precisely ^ Thus, while 

definition I gives no suggestion of hyperuniform behavior, 
definition II looks convincingly hyperuniform out to rel¬ 
atively large length scales R. The dramatic difference in 
var(0/j) between the two dehnitions of the local packing 




FIG. 11. (color online) R^var{(j)ii) vs window radins R for 
systems with N = 8192 to 65536 particles, at stress (a) p = 
0.0001373 and (b) p = 0.0183105. The power law decay at 
small R is indicated. Definition I of Eq. Q for the local 
packing fraction (j){r) is used. 


fraction (p{r) is quite puzzling given the complete agree¬ 
ment of the corresponding x(q) for the two definitions at 
small |q|, as seen in Fig. We can explain the reason 
for this difference in behavior as follows. 



10 ^ 100 


10 ^ 100 


FIG. 12. (color online) R?M&r{(j)R) vs window radius R for 
systems with N = 8192 to 65536 particles, at stress (a) p = 
0.0001373 and (b) p = 0.0183105. The power law decay at 
small R is indicated. Definition II of Eq. 0 for the local 
packing fraction (f){r) is used. 


From Eqs. (14) and ( |19[ ) 
tween var((^/j) and x(qjas, 


we can write the relation be- 


var((^fl) = ^ XI x(q)/^(|q|-R), (27) 

q#0 


with f{y) as defined in Eq. ( [T4| . Assuming that x(q) 
depends only on |q| due to the average isotropy of the 
system |31], we can integrate over the direction of q to 
get for our two dimensional system, 

Yav{(j)R) = ^ X (28) 

y q^O 


with q = 27TnjLy. In Fig. 13 we plot f^{y) vs y on a 
log-log scale. We see that forTarge y, it oscillates with a 
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1/y^ envelope. For an infinite system, if y(q —>■ 0) is a 
finite constant, then at sufficiently large R a dimensional 
analysis implies that var(^fl) must scale as How¬ 

ever for finite R, and in finite systems where the sum 
on q is discrete, the behavior of var( 0 /j) can depend in 
detail on the behavior of x(q) at large q; if the sum in 
Eq. (28) is dominated by the large q terms, then we may 
find va,T{(j)ji) ~ 1/i?^ because of the l/{qR)^ dependence 
of f^{qR) and not because of any hyperuniformity of the 
system. 



FIG. 13. (color online) Plot of f^{y) vs y, with the function 
f{y) as defined in Eq. ( |14[ ). The envelope of the oscillations 
at large y decays as 1/yT 



FIG. 14. (color online) Packing fraction fluctuation xilf) 
computed using definition I of Eq. ([^, and definition II of 
Eq. (0, as well as the structure function S{qy), for a wide 
range of wavevectors up to g = 10. Results are shown for 
a system oi N = 32768 particles at our lowest stress p = 
0.0001373. 


that the reason for the dramatic difference in 
comparing definition I with definition II, is the influence 
of fluctuations at moderately large q, which persist even 
to large R. We conclude that x(q), rather than var(^fl;), 
is the better measure to use to check for hyperuniformity 
in our two dimensional system. 


The behavior of var((/)j{) on observed length scales R 
can thus be determined by the behavior of x(ct) s-f large 
wavevectors q with |q| > ir/R. In Fig. 14 we plot x('?y) 
for both definition I and definition II, as well as the struc¬ 
ture function S{q'y), for a much wider range of wavevec¬ 
tors, 0 < 9 < 10, than in previous plots. We show results 
for N = 32768 at our lowest stress p = 0.0001373. As 
before, we see that x(<zy) foi' the two definitions agree 
perfectly at small q, but then separate when g ^ 1 . 
Moreover, xiQV) tor definition II becomes roughly equal 
to S{qy), and over two orders of magnitude larger than 
that for definition I, when g ^ 5. Thus the contribu¬ 
tion to var( 0 fl') from x( 9 y) at large g should be expected 
to be more significant for definition II as compared to 
definition I. 

To check this, we compute var((^ij) by explicitly sum¬ 
ming the series in Eq. (28), using the data for x('Zy) 
from Fig. 14 Our results for R^vaT{<pji) vs R are shown 


in Fig. |15| for N = 32768 at p = 0.0001373. We com¬ 
pare these results to the direct computation of var( 07 j) 
as shown in Figs. [10, and[T0.. For definition I we And ex¬ 
cellent agreement between the series and the direct com¬ 
putation when we sum the series up to g = 10. For defi¬ 
nition II we And that we must sum even more terms, up 
to g = 50, in order to get reasonable agreement. We thus 
see again that the large g (i.e. small length scale) fluctu¬ 
ations are larger, and so contribute more to var(^ij), for 
definition II than for definition I. Our results in Fig. |0| 
show that the computation of var(0fl') in Fig.[0]is indeed 
consistent with our computation of xiQy) in Eig. 0 and 



FIG. 15. (color online) R^var{(j)ii) vs R, comparing the data 
of Figs. and |12^ (direct computation) with the result from 
summing xiiy) <1 in Eq. ( |28| ), for both definition I and 
definition II of the local packing fraction Results are 

for a system with N = 32768 particles at our smallest stress 
p = 0.0001373. 


V. DISCUSSION AND CONCLUSIONS 

In this work we have considered the fluctuations of den¬ 
sity and packing fraction in mechanically stable packings 
of bidisperse frictionless particles at finite stress. A dis¬ 
tinguishing feature of our work is that we simulate at 
fixed isotropic global stress, rather than at fixed packing 
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fraction. We investigate states above the jamming tran¬ 
sition, in contrast to earlier works [aHii] that considered 
the case of packings exactly at the jamming (j}j. 


A. Comparison to previous works 

Berthier et al. [S] considered both an experimental 
two dimensional system of IV = 8000 bidisperse particles, 
and a numerical three dimensional system of = 64000 
soft-core particles of varying size dispersities. This cor¬ 
responds to system lengths of roughly 90 and 40 parti¬ 
cle diameters in the experimental and numerical systems 
respectively. The experimental system was jammed un¬ 
der slow compression. The numerical configurations were 
created starting from jammed states above (f>j, and then 
slowly decompressing until the system unjammed, fol¬ 
lowed by a slower recompression until the system jammed 
again, as measured by a finite energy per particle of 
order 10“^^. In both experimental and numerical sys¬ 
tems, a xt(<?) is observed that appears to linearly de¬ 
crease towards zero. However, in both cases the mea¬ 
sured xt{q) only extends down to roughly qd « 0.15, 
where d is the average diameter of the particles, and the 
data is quite scattered below qd « 0.5. These results thus 
give evidence for hyperuniformity out to length scales 
£/d ^ 27r/0.5 « 12, but not necessarily on longer length 
scales. 

Zachary et al. mm considered a much larger nu¬ 
merical system in two dimensions, with up to iV = 10® 
particles (and so system length of roughly 1000 particle 
diameters). They used a bidisperse system with particles 
of different shapes, with size ratio 1.4 and a concentration 
of small particles Xs = 0.75, and large particles Xb = 0.25. 
They used the Lubachevsky-Stillinger algorithm [33] to 
generate their jammed states. This is an event-driven 
molecular dynamics for elastic hard-core particles, where 
particles are inflated at a prescribed rate from an ini¬ 
tial thermally equilibrated dilute state so as to rapidly 
quench the hard-core gas into a thermal glassy state. 
The particle inflation continues until the system seizes 
up into a strictly jammed state that they denote as max¬ 
imally random jammed (MRJ). They measure both x{d) 
and var((/)fl) (using definition I) in the MRJ state, and for 
circular particles they find strong evidence from var( 07 j) 
for hyperuniformity out to length scales 2R/d ^ 80. 

More recently, other works have reconsidered hyper¬ 
uniformity in systems of monodisperse spheres in three 
dimensions, and have considered behavior approaching, 
rather than strictly at, (j)j. Recall, for monodisperse sys¬ 
tems, hyperuniformity is indicated by the behavior of the 
structure function, S'(q) ~ |q| as |q| —)■ 0. Hopkins et al. 

using the same protocol as Zachary et al. dnuni for 
a system with N = 10® particles (system length roughly 
100 particle diameters), measure 5'(q) at various (j) ap¬ 
proaching (pj from below. They find 5'(q) ~ aq-\-b, with 
6 —>■ 0 as (/) —>■ and from this extract a length scale 

^ ^ that diverges as jamming is approached and 


the system becomes hyperuniform. 

Ikeda and Berthier [Mj study N = 512000 monodis¬ 
perse soft-core particles in three dimensions. Starting 
from a random configuration of particles in a fixed cubic 
box at packing fraction (j) = 0.8, well above (j)J ~ 0.646, 
they use the FIRE algorithm |5^ to minimize the interac¬ 
tion energy and obtain a mechanically stable state. They 
then decrease the particle density in small steps, energy 
minimizing at each step, to obtain configurations span¬ 
ning a range of packing fractions from (j) = 0.8 to just 
above the jamming (j)j. Their results are averaged over 
8 independent starting configurations. Computing S'(q) 
they find, for all but their largest value of 4> = 0 . 8 , that 
data at the different (j) essentially overlap and are linear 
in q, as expected for a hyperuniform system, over an ex¬ 
tended range of 0.4 < q <7. However at their smallest 9 , 
they find 5'(q) saturates to a finite value ~ 10 “®, similar 
in magnitude to what we have found in the present work 
for x(q); they, however, see a plateau in S'(q) at small 
q rather than the minimum that we find in x(q)- Ikeda 
and Berthier thus conclude that hyperuniformity is only 
weakly dependent on packing fraction </>, but persists out 
to only a finite length scale « ISd. Ikeda and Berthier 
further find that this behavior is stable to the addition 
of small finite thermal fluctuations. 

The above results, combined with our own, suggest 
that mechanically stable jammed packings above (jj do 
not display hyperuniform fluctuations of the packing frac¬ 
tion out to arbitrarily large length scales, but are hyper¬ 
uniform only out to a finite £*{(j)) that is weakly depen¬ 
dent on (j) and does not appear to diverge as 4> ^ (pj from 
above. However the results of Zachary et al. [Mini and 
Hopkins et al. [T 2 | suggest that hyperuniformity may ex¬ 
ist in hard-core particle systems, when compressed to pj 
from below. It may therefore be that the presence or ab¬ 
sence of hyperuniformity out to arbitrarily large length 
scales depends on the specific protocol used to construct 
the jammed state at pj. We also cannot rule out the pos¬ 
sibility that hyperuniformity may still exist in jammed 
packings above pj, but restricted to a region closer to 
pj than we have been able to explore in this work. 


B. Alternative ensembles 


To check how sensitive our results for p > pj are to the 
particular system we have used above, we have consid¬ 
ered two other ensembles. The first is to use a Hertzian 


interaction, with a = 5/2 in Eq. (22 1 , in place of the har¬ 


monic interaction. All other details of the system remain 
the same. In Eig. [^we plot x(<zy): computed accord¬ 
ing to definition I of Eq. ([^, vs g for the four lowest 
p that we used for the harmonic interaction. We use a 
system size with N = 32768 particles. Eor the Hertzian 
interaction, pressure is expected [ 2 | to scale with packing 
fraction according to p ^ (p — pjY^"^, hence the (p) for 


the Hertzian system (shown as the inset to Fig. 16) is 


larger than that of the harmonic system at equal values 
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FIG. 16. (color online) Fluctuation in packing fraction xiw) 
vs q, using the definition I of Eq. Q for the local packing 
fraction (()(r), for the case of a Hertzian interaction (a = 5/2 
in Eq. (22l). Results are shown for several different values 
of the stress per particle p, for a system with N = 32768 
particles. The inset shows the average packing fraction (</) as 
a function of p. 


of p. The Hertzian interaction further differs from the 
hamonic in that for the Hertzian system the bulk mod¬ 
ulus vanishes continuously as ^ ^ t/fj from above, while 
for the harmonic system the bulk modulus approaches a 
finite constant as (p ^ (j)J from above, and then jumps 
discontinuously to zero below pj [5]. Nevertheless, we 
find that xi^y) for fh® Hertzian system is qualitatively 
the same as for the harmonic case, with a well defined 
minimum that does not appear to be moving to smaller 
q as p decreases. 

The second is to consider the harmonic interaction, 
but to obtain our configurations by quenching at fixed 
(j) within a fixed square box. Such constant cp ensembles 
have usually been used in earlier works [SiiaiSl]. Un¬ 
like the constant stress ensemble, where configurations all 
have the same p and so can be viewed as all at the same 
distance from the jamming transition p — 0, the constant 
p ensemble has a fluctuating p and so different configura¬ 
tions i are at different distances from their configuration 
specific jamming transition pji UM- The constant p 
ensemble in a fixed box also allows there to be a finite 
residual shear stress in the quenched configuration |18j . 
In Fig. 17 we plot x(gy), computed according to defini¬ 
tion I of Eq. ([^, vs g for the case p = 0.8422 close to 
pj « 0.84159. We use a system size with N = 32768 
particles. Again we see qualitatively the same behavior 
as before. 

Note, the fixed value oi p = 0.8422 in Fig. was 
chosen as it is equal to the (p) for a system with p = 
0.0002747 in the fixed stress ensemble (our next to lowest 
value of p). However in the fixed p = 0.8422 ensemble, we 
find that the average stress per particle is (p) = 0.000180, 
lower than the corresponding value in the fixed stress en¬ 
semble. This suggests that the jamming density pj of the 
constant p ensemble is slightly larger than the jamming 


FIG. 17. (color online) Fluctuation in packing fraction x(gy) 
vs q, using the definition I of Eq. <© for the local packing 
fraction p{r), for the case of the harmonic interaction in an 
ensemble at fixed global packing fraction p = 0.8422. Results 
are shown for a system with N = 32768 particles. 


density of the constant stress ensemble. That is consis¬ 
tent with our estimate of pj « 0.84159 for the constant 
stress ensemble from Fig. as compared with the esti¬ 
mate of pj « 0.84177 for the constant p ensemble from 
Ref. [H]. We also note that the width of the distribu¬ 
tion of p found in this fixed p ensemble is rather large, 
V ~ {P)'^/P = 0.40, while the corresponding width 
of the distribution of the r esidu al deviatoric stress per 
particle a is rather small, \fW^lp = 0.00053. 


C. Rattlers and polydispersity 

It has been suggested giMi that rattlers may play a 
role in the breaking of hyperuniformity on large length 
scales. Rattlers result when a particle has an insufficient 
number of contacts to constrain its motion in all direc¬ 
tions. Determining the number of rattlers according to 
the method described in Sec. |IV A| in Fig. [T^ we plot 
the fraction of particles that are rattlers (Afrattiers)/-^ vs 
the stress per particle p. For the harmonic interaction, 
we plot results for systems with N = 8192 to 65536 par¬ 
ticles. For the Hertzian interaction, we plot results for 
N = 32768 only. We see that (A^rattiers)/-^ is indepen¬ 
dent of the system size N^ and decreases with increasing 
p. For the harmonic interaction, (Ni.attiers)/A^ changes 
by an order of magnitude over the range of p we study. 
If rattlers were responsible for the breaking of hyperuni¬ 
formity, we might expect that the length I* to which 
hyperuniformity extends should increase as the density 
of rattlers decreases, i.e. as p increases. However our 
results in Fig. show exactly the opposite trend; the q* 
that locates the minimum of xi^v) increases slightly with 
increasing p, and so i* = ‘2 ,'k jq* decreases with increasing 
p. Our results thus provide no obvious relation between 
rattlers and the breaking of hyperuniformity. 

Another possibility that might lead to the breaking of 
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FIG. 18. (color online) Fraction of the particles that are rat¬ 
tlers, (Afrattiers)/^'^ VS stress per particle p = Fn/N. For 
the harmonic interaction we show results for systems with 
N = 8192 to 65536 particles; the points for different N in 
the figure overlap each other. For the Hertzian interaction we 
show results for N = 32768. 


hyperuniformity is suggested |36j by the work of Dreyfus 
et al. m- Their work is primarily concerned with the 
detection of hyperuniformity in experimental systems, 
where particles are polydisperse, and the exact size of 
individual particles is not a priori known but must be 
determined by optical measurements. Errors in the de¬ 
termination of the exact particle sizes were found to re¬ 
sult in an apparent breaking of hyper uniformity at small 
wavevectors (large length scales). As an extreme example 
of this effect, one can consider the error introduced if, in 
a bidisperse or polydisperse system, one approximated all 
particles as having the same average size. In that approx¬ 
imation, the packing fraction fluctuation x(q) J^st be¬ 
comes proportional to the structure function S (q) , which 
clearly does not show hyperuniformity at small q, as seen 
in Fig. 

In our bidisperse particle simulations, we of course 
know the position and size of each and every particle 
exactly. Nevertheless, an effective polydispersity may be 
viewed to arise from the following effect. Both our defi¬ 
nitions I and II count each particle with a weight equal 
to the area of the particle in isolation. However in our 
jammed packings, particles in contact necessarily have 
some amount of overlap. Our definitions I and II there¬ 
fore count this overlap area twice, once for each particle. 
One might imagine that a more “correct” definition of the 
local packing fraction should count this overlap area only 
once, dividing it proportionally between the two contact¬ 
ing particles. For example, as sketched in the inset to 
Fig. particle i should have a weight equal to only 
the shaded area, rather than the full area of the corre¬ 
sponding circle. If Vi = is the area of the circle 

of particle f, then the weight with which particle i en¬ 
ters the local packing fraction should instead be taken as 
Vi = Svij, with dvij the area subtracted due to the 

overlap with particle j. The weights Vi are therefore poly¬ 
disperse, depending on the varying overlaps in the sys¬ 


tem. If one computes x(q) using the bidisperse weights 
Vi rather than the more correct polydisperse weights Vi, 
it could lead to a breaking of hyperuniformity that is only 
apparent, i.e. a consequence of using incorrect weights. 

However, if dij = {di + dj)/2 — Vij is the overlap length 
of the contact, then Svij /vi (x {6ij/di)^^^ (X where 
the last result follows since the pressure p ~ {Sij) for the 
harmonic interaction potential. Thus this effect should 
vary with the pressure and vanish continuously as p —)■ 0, 
as one approaches the jamming transition. To test this 
notion, we have therefore computed xidY) according to 
definition H of Eq. 0, but using the weights Vi as de¬ 
scribed above, computed exactly for each particle ac¬ 
cording to its own specific overlaps. We use definition 
H since it is easier to implement than definition I, in the 
case where each particle has a unique, nonsymmetric (i.e. 
circle minus overlaps), shape. However we expect from 
Fig. that xidy) be identical for definitions I and H 
at the small q of interest. In Fig. 19 we plot the resulting 
x{qy) for our largest system with N = 65536 particles, 
at both our smallest and largest values of p. We compare 
the xidy) obtained from using the original weights Vi (de¬ 
noted as “counting overlaps twice”) with that using the 
new weights Vi (denoted as “counting overlaps once”). 

At the largest p = 0.0183105, we see a clear shift be¬ 
tween the results from the different sets of weights, how¬ 
ever the qualitative behavior remains the same, with a 
clear minimum at the same q*, and x increasing as q 
decreases below q*. At our smallest p = 0.0001373, how¬ 
ever, the results from the two sets of weights are essen¬ 
tially equal. Thus taking overlaps into account does not 
result in a restoration of hyperuniformity on large length 
scales, and the insensitivity of our results to the differ¬ 
ent choices of weights at our smallest p is yet another 
indication that our smallest pressures are, by all relevant 
measures, quite close to jamming. 

Finally, it is interesting to note that the experiments 
on PINIPAM microgel particles, reported on in Dreyfus 
et al. may actually correspond more closely to our 
conclusions than to the claim in favor of hyperuniformity. 
As these authors note, the small q behavior of xid) for 
PINIPAM, shown in the inset to their Fig. 9b, does not 
suggest hyperuniformity; indeed it is qualitatively sim¬ 
ilar to what we see in our Fig. 7. However the q* at 
which xid) bas its minimum is so much larger in the ex¬ 
periments of Dreyfus et al. than what we find here, that 
in their case it may well be an artifact of system size, 
such as Dreyfus et al. claim. However, if we consider the 


real space decay of with R, our results in Fig. II 


for definition I (corresponding to the usage in Dreyfus et 
al.) show that the initial decay, before the minimum is 
reached, is va,r{(j)]i) ~ R~^, with A « 2.4 for our smaller 
p, and A « 2.3 for our larger p. This is not far from the 
value A « 2.2 reported in Dreyfus et al. for a similar 
range of R, using their j-PSR reconstruction as shown 
in their Fig. 9a. Yet in our case, our A > 2 does not 
demonstrate that the system is hyperuniform; we see hy¬ 
peruniformity is broken only by looking at larger length 
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FIG. 19. (color online) Packing fraction fluctuation xiiy) ''S 
q, using definition II of Eq. 0 for a system of = 65536 par¬ 
ticles, at two different values of the stress per particle p. We 
compare results where the weight of each particle is taken as 
the circular area of the isolated particle (denoted as “count¬ 
ing overlaps twice”), vs where the weight of each particle is 
taken as the non-overlapping part of that circular area, as 
illustrated by the shaded region for particle i in the inset (de¬ 
noted as “counting overlaps once”). A small difference is seen 
between these two sets of weights at the larger value of p, but 
not at the smaller value. 


scales. This comparison thus suggests that the PINIPAM 
experiments may actually be above the jamming (j)j, and 
are not inconsistent with the absence of hyperuniformity 
on long length scales. 
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APPENDIX 


In this appendix we provide some further details about 
the minimization procedure of Sec. mil that we use to 
obtain our mechanically stable configurations at fixed 
isotropic stress. 

Since our minimization procedure is carried out at 
fixed total system stress Eq,/ 3 , the system box parameters 
Lx, Ly and 7 (see Fig. will vary from specific mini¬ 
mized configuration to configuration. In Fig. |20| we show 
the extent of these variations for the different system sizes 
N = 8192, 16384, 32768 and 65536. In Fig [2^ we show 
the relative fluctuations in box lengths, ^yvai(Lx)/{Lx) 


and y/var{Ly)/(Ly) vs the stress per particle p = Tn/N. 
Solid symbols are for Lx while open symbols are for Ly. 
Since the system is on average isotropic, we expect the 
fluctuations in Lx and Ly to be equal, and we indeed find 
that to be so. The fluctuations are also found to scale 
as 1/Vn, as would naively be expected. In Fig. 


20 3 we 


show the fluctuations in the dimensionless skew parame¬ 
ter, .^var( 7 ) vs p. The size of the fluctuations in 7 are 
slightly larger but comparable to the fluctuations in the 
box lengths. Again we find that the fluctuations scale as 
l/y/N. We also note that, as expected, the average skew 
( 7 ) = 0 within the estimated statistical error, as shown 
in Fig. [2T| 
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FIG. 20. (color online) (a) Relative fluctuations in the sys¬ 
tem box lengths Lx and Ly vs stress per particle p, for system 
sizes N = 8192 to 65536. Solid symbols show the fluctuations 
in Lx, while open symbols show Ly. (b) Fluctuations in the 
dimensionless box skew parameter 7 vs p. See Fig. [^for the 
definition of parameters Lx, Ly, 7 . In both cases the fluctu¬ 
ations scale as 1/%/lV. 



10 -“ 10-3 n 10 '^ 


FIG. 21. (color online) Average box skew ( 7 ) vs stress per 
particle p, for system sizes N = 8192 to 65536. Error bars 
represent one standard deviation of estimated statistical error, 
showing that ( 7 ) = 0 within the estimated errors. 


Our minimization procedure necessarily produces the 
desired isotropic stress configurations only to a certain 
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numerical accuracy. We now provide details of the degree 
of that accuracy. 

We first look at how well our procedure produces a 
packing with the desired isotropic global stress tensor, 
Sa /3 = r N^ap. We compute the global stress tensor Tiap 
for our minimized configurations using the usual formula 
[5] for a static frictionless system, 


we plot the average magnitude of the net force, nor¬ 
malized by the average magnitude of the contact force, 
(|Fi|)/(|Fy|), vs the stress per particle p = Tf^/N, for 
system sizes N = 8102 to 65536. We see that the resid¬ 
ual net force on a particle at the end of our minimization 
procedure is less than 0.05% of the average contact force. 
This decreases as either p or N increases. 


^a/3 — ^ ^ jQ ^ijp j 

i<3 

where Yij = Yj — is the center-to-center displacement 
from particle i to particle j, F^- = —dVijjdYi is the con¬ 
tact force on i due to j, and the sum is over all distinct 
pairs of particles in contact. We then define three mea¬ 
sures of the deviation of our minimized stress from the 
isotropic target value. 


5i = 
62 = 


2{^xx + ^yy) ~ TAf]^) 


Tat 


63 = 


i N 

_ 


■ N 


(30) 

(31) 

(32) 


i5i measures the relative spread in the trace of "Sap about 
the target value S 2 measures the relative spread in 
anisotropy of the diagonal elements of ; and ^3 mea¬ 
sures the relative spread in the off-diagonal elements of 
Tiap. In Fig. 22 we show our results for di, 62 and 63 vs 


p = Tn/N, for systems sizes N = 8192 to 65536. We see 
that (5i is less than 0 . 01 %, while S 2 and ^3 are less than 
0.004%, indicating a high accuracy in the desired stress 
tensor. In all cases the accuracy improves as the stress 
per particle p increases, and as the number of particles 
N increases. 



FIG. 23. (color online) Average magnitude of the residual net 
force on a particle F^, normalized by the average magnitude of 
the contact force F^, at the termination of our minimization 
procedure. Results are plotted vs the stress per particle p = 
r n/N for systems sizes N = 8192 to 65536. 


Fig. |23| showed the average net residual force on parti¬ 
cles. In Fig. we show the distribution of such forces, 
7^(|Fi|/(|Fij|)) vs |Fi|/(|Fij|), for different system sizes 
N = 8102 to 65536, at our (a) smallest p = 0.0001373 and 
(b) largest p = 0.0183105. We see that the large force tail 
grows as N increases, but shrinks as p increases. For our 
largest system, N = 65536, at our lowest stress per parti¬ 
cle, p = 0.0001373, there exist a very few particles whose 
net force is comparable to the average contact force. 




FIG. 22. (color onli ne) Ac cur acy pa ram eters (a) 5i, (b) S 2 , 
and (c) 53 of Eqs. (301, (311 and (32 1 , that measure the 


relative deviations of the stress tensor from the target 
isotropic F N^ap, vs stress per particle p = F n/N, for system 
sizes N = 8102 to 65536. 


Next we look at how well our procedure produces a me¬ 
chanically stable packing in which the net force on each 
particle vanishes. The net force F^ on particle i is just 
the sum over its contact forces 



- Z) 7 ■ 


In Fig. 23 


FIG. 24. (color online) Distribution 'P(|Fi|/(|Fijj)) of the 
residual net force on particles |Fi|, normalized by the aver¬ 
age magnitude of the contact force |Fijj, vs |Fi|/(|Fijj), for 
different system sizes N = 8102 to 65536, at our (a) smallest 
p — 0.0001373 and (b) largest p = 0.0183105. 


The average residual force (|Fi|), and the large force 
tail of the distribution, is controlled by the accuracy 
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parameter e that determines when we stop our mini¬ 
mization procedure, (Ui — C/i+ 5 o)/C/i+ 5 o < e. In the 
body of this work, and i n th e above results, we have 
used £ = 10“^°. In Fig. 25 we show the distribution 
■PdFil/dFijD) for several different values of the accu¬ 
racy parameter 10 ”^° < e < 10 “^, for our biggest sys¬ 
tem N = 65536 at our lowest stress p = 0.0001373. We 
see that as £ decreases, the average net force and the 
large force tail decrease. Thus, as would be expected, 
decreasing e improves the accuracy of force balance on 
the particles in our minimized configurations. 





FIG. 25. (color online) Distribution PdFil/dFij|)) of the 
residual net force on particles |Fi|, normalized by the average 
magnitude of the contact force |Fij|, vs |Fd/{|Fij|), for dif¬ 
ferent values of the accuracy parameter e that terminates onr 
minimization of U. Results are for a system of size N = 65536 
at p = 0.0001373. Data labeled “10“^°*” are results for a 
minimization of U to accuracy e = 10~^'^, followed by an 
additional minimization of interaction energy U to accuracy 
10 “^° while holding the box parameters constant. 


We have attempted to improve upon the accuracy of 
force balance by adding a separate step of minimization 
in which, after the above criterion on U is met, we then 
hold the box parameters Ly and 7 constant while 
adjusting the particle positions to minimize the interac¬ 
tion energy, ([/* - I/^+soj/t/i+so < 10“^°. The resulting 
distrib utio n of net residual forces on particles is shown 
in Fig. [ 2 ^ labeled as “10“^°*”. We find a significant re¬ 
duction in the net force, with the average dFd)/dFy|) 
decreasing roughly by a factor of 100. However we also 
find that the accuracy of the system to have the desired 
target global stress decreases, with the parameters 5i of 
Eqs. (30||^l increasing roughly by a factor 10. We have 


not tried to optimize the sequence of minimizing U and 
U as we have found our results for x(q) to be insensitive 
to this additional step of minimization (see below), and 
so we have not used it for the results presented elsewhere 
in this paper. 


Finally, to determine whether the accuracy parameter 
£ = 10 “^° used in this work is sufficient for our needs, 
we now check the sensitivity of x(gy) to the value of 
£. In Fig. we plot x(<?y) vs q (using definition I of 
Eq. for ^(r)) for different values of £ = 10“® to 10“^°, 
for a system with N = 65536 particles (we consider our 
largest system since that has the force distribution with 
the largest tail at large |Fi|). We show results for our 
smallest and largest values of the stress per particle p = 
Tff/N. We see that if £ is too large, the results at small q 
are clearly dependent on e. But as £ decreases, our results 
converge to a fixed £-independent curve. For the smallest 
p = 0.0001373 this happens for e < 10“®, while for our 
largest p = 0.0183105 we have convergence for £ < 10“®. 
For the lowest p in panel (a) we also show results for 
the case where £ = 10 “^° and we add the second step of 
minimization described above, in which we fix the box 
parameters and only move particle positions to minimize 
U. This data is labeled as “10~^°*” in the figure. We see 
that this additional step of minimization does not result 
in any noticeable change in xi^y)- We thus conclude 
that using e = 10 “^° with a single step minimization of 
U gives sufficient accuracy for our needs. 
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(b) our largest value p = 0.0183105. We compare the values of 
x(gy) obtained when using different values of the parameter 
e that determines the stopping criterion for our minimization 
of U. Data labeled “10“^°*” in panel (a) are results for a 
minimization of U to accuracy e = 10“^°, followed by an 
additional minimization of interaction energy U to accuracy 
10~^'^ while holding the box parameters constant. 
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